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Abstract. We investigate the cosmological constraints on axion models where the domain 
" wall number is greater than one. In these models, multiple domain walls attached to strings 

are formed, and they survive for a long time. Their annihilation occurs due to the effects of 
explicit symmetry breaking term which might be raised by Planck-scale physics. We perform 
three-dimensional lattice simulations and compute the spectra of axions and gravitational 
waves produced by long-lived domain walls. Using the numerical results, we estimated relic 
density of axions and gravitational waves. We find that the existence of long-lived domain 
walls leads to the overproduction of cold dark matter axions, while the density of gravita- 
tional waves is too small to observe at the present time. Combining the results with other 
observational constraints, we find that the whole parameter region of models are excluded 
unless an unacceptable fine-tuning exists. 
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1 Introduction 

Axion [1] is a hypothetical particle which arises as a consequence of the Peccei-Quinn (PQ) 
mechanism [2], the most attractive solution to the strong CP problem of quantum chromo- 
dynamics (QCD) [3, 4]. It arises as a Goldstone boson when a global C/(1)pq symmetry (so 
called PQ symmetry) is spontaneously broken at some high-energy scale denoted as F a . The 
phenomenological searches indicate that axions are "invisible" due to the weak coupling with 
matter which is suppressed by 1/F a . This invisibleness is realized by two kind of models. One 
is the KSVZ models [5], which is realized by introducing additional quarks to the standard 
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model. Another is the DFSZ models [6], which is realized without introducing additional 
fermions by assuming two Higgs doublets. These invisible axions lead to rich implications to 
astrophysics and cosmology, such that they are good candidates of dark matter filled in the 
universe [7]. 

The remarkable feature of the axion models is that they predict the formation of topo- 
logical defects in the early universe [see [8] for reviews]. When the PQ symmetry is spon- 
taneously broken, vortex-like defects, called strings, are formed. These strings are attached 
by surface-like defects, called domain walls, when the temperature of the universe becomes 
comparable to the QCD scale. The structure of the domain walls is determined by an integer 
number Apjw which is referred as the "domain wall number" . The value of iVow is related 
to the U (1)pq-SU (3) C -SU (3) c anomaly which is given by the formula [4, 9, 10] 



(1.1) 



2 £ TrQ PQ fe)T a 2 ( % ) - 2 TrQ PQ ( qi )T^ qi 

i=L i=R 

where Qpq(%) is the J7(1)pq charge for quark species qi, Yli=L{R) represents the sum over 
the left(right)-handed fermions, and T a are the generators of SU(3) C normalized such that 
TiTaTf, = I5 a b where 1 = 1/2 for fundamental representation of SU(3) C . Here we choose the 
unit Qpq = 1 for minimum non-vanishing magnitude of PQ quantum numbers. In the KSVZ 
models we obtain Adw — 1) while in the DFSZ models it is predicted that Adw = 27V fl 
where N g is the number of generations. As we discuss in Sec. 2, the cosmological scenario is 
different between models with Afow = 1 and Afow > 1- Recently, we studied the cosmological 
implication of the models with Afow = 1 through the careful numerical simulations [11]. In 
this paper, we turn our research into the models with A^dw > 1- 

For the case A^ow > 1 , it is well known that domain walls are stable and they eventually 
overdose the universe, which conflicts with the standard cosmology [9, 12]. So far, several 
solutions to this problem are proposed. One is to introduce a small symmetry breaking 
parameter, called the bias [13], in the Lagrangian, which makes the walls unstable and leads 
to their annihilation at late times. Another solution is to embed the discrete subgroup 
of J7(1)pq into another continuous group [14]. In such models, the degenerate vacua are 
connected by symmetry transformations, leading to the same cosmological property with 
the case for Afow = 1- It is also possible to avoid the domain wall problem by assuming 
that inflation has occurred after the PQ symmetry breaking. In this case, the inflationary 
scale is constrained by the observation of isocurvature fluctuations in the cosmic microwave 
background (CMB) [15]. Among these possibilities, here we focus on the first solution, the 
biased domain walls. 

It was pointed out that the bias term, which causes the annihilation of domain walls, 
can be found in effects of gravity [16]. If we regard the low energy field theory as an effective 
theory induced by a Planck-scale physics, we expect that the global symmetry is violated 
by higher dimensional operators suppressed by the Planck mass Mp. This small symmetry 
violating operators cause the decay of domain walls. However, in the case of PQ symmetry, 
the situation is more complicated. It was argued that these Planck-scale induced terms easily 
violate the PQ symmetry and cause the CP violation [17-19]. In particular, the operators 
with dimension smaller than d = 10 are forbidden [18] by a requirement of CP conservation. 
If it is true, the biased domain walls should be long-lived, and disappear at late time due to 
the effect of highly suppressed operators. 

In our previous numerical study [20] , we found that the late-time annihilation of domain 
walls actually occurs. We also pointed out that the gravitational waves produced by domain 
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walls may be detectable in the future experimental searches. However, the previous results 
strongly depended on the production rate of axions and gravitational waves, which was 
treated as free parameters. In this work, extending the previous study [20], we aim to 
determine the amount of radiated axions and gravitational waves precisely. This requires the 
estimation of the spectra of radiated particles from domain walls, which can be investigated 
by performing field-theoretic numerical simulation of scalar fields. 

In this paper, we study the formation and evolution of string-wall networks in the 
models with Adw > 1 an d calculate the spectra of axions and gravitational waves radiated 
from them. This work might be regarded as the first study in which the particle radiations 
from long-lived domain wall networks are estimated quantitatively. Although we focus on the 
axionic models, part of our findings might be applicable to other particle physics models which 
predict the existence of long-lived domain walls, such as domain walls formed by discrete R- 
symmetry breaking [21], thermal inflationary models [22], and Z% symmetry breaking in 
NMSSM [23]. 

The organization of this paper is as follows. We introduce the model and review the 
property of topological defects in Sec. 2. In Sec. 3, we show the results of numerical simula- 
tions. Using the numerical results, we calculate the relic density of axions and gravitational 
waves in Sec. 4. We show that the overabundance of axions produced by domain walls give 
severe constraints on the model parameters. We also review other observational constraints. 
Finally, we make a summary in Sec. 5. The analysis method that we used in the numeri- 
cal studies are described in Appendices A, B, C, and D. We make a brief comment on the 
determination of the surface mass density of domain walls in Appendix E. 

In this paper, we work in spatially flat Friedmann-Robertson- Walker (FRW) universe 
with a metric 

ds 2 = dt 2 - R 2 (t)[dx 2 + dy 2 + dz 2 ], 

where R(t) is the scale factor of the universe. We denote the cosmic time as t and the 
conformal time as r, where dr = dt/R(t). A dot represents a derivative with respect to the 
cosmic time, while a prime represents a derivative with respect to the conformal time, i.e. 
-=d/dt, and ' = d/dr. 

2 Axion cosmology and topological defects 
2.1 Model of domain walls bounded by strings 

We consider the model of a complex scalar field $ with the Lagrangian density given by 

£ = \d»w* - v($), (2.i) 

where the scalar potential is given by 

V(*) = j(m 2 - V 2 ) 2 + - cosiW?)- (2.2) 

4 JV DW 

Here, 9 is the angular part of the scalar field (i.e. = l^le 40 ). The first term in Eq. (2.2) 
is the usual mexican-hat potential, which causes spontaneous breaking of C/(1)pq symmetry. 
In the early universe, the C/(1)pq symmetry is broken when the temperature of the universe 
falls below the energy scale ~ r/. At this time, the scalar field gets a vacuum expectation 
value K^)] 2 = T] 2 , and cosmic strings are formed. Then, the angular part of the scalar field 
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can be identified as axion field a, <1> = ?7e m / r? . The second term in Eq. (2.2) may arise due 
to the non-perturbative effect of QCD. Because of the presence of this term, the C/(1)pq 
symmetry is explicitly broken into its subgroup Zn dw . This term becomes non-negligible 
when the mass of the axion m becomes greater than the friction term which comes from the 
cosmic expansion, m > H, where H is the Hubble parameter. Then, the axion field begins 
to roll down to one of Adw degenerate minima represented by 

Ink 

fc = 0,l,...,iVDw-l- (2.3) 



At this stage, domain walls are formed, separating the Afow vacua [9]. These domain walls 
are attached to cosmic strings around their boundaries [8]. 

We note that there are two different energy scales in this model. One is the scale of the 
PQ symmetry breaking, given by the axion decay constant 

F a = 11 (2.4) 

Phenomenologically, it takes a value F a ~ 10 9 -10 12 GeV (see e.g. [3]). Another is the QCD 
scale, Aqcd — C>(100)MeV, which induces the second term of the effective potential (2.2). 
The mass of the axion is given by m ~ O(0.1)AQ CD /F a . Since there is a large hierarchy 
between these two scales F a 3> Aqcdj the height of the potential barrier of the first term 
~ r/ 4 in eq. (2.2) is much larger than that of the second term ~ Aq CD . 

Because of the existence of this large hierarchy, the formation of strings occurs much 
earlier than the formation of domain walls. Once strings are formed, they evolve into the 
"scaling" regime, in which the population of strings in the Hubble volume tends to remain in 
0(1). In order to keep this scaling property, the strings lose their energy by emitting closed 
loops, and these loops decay by radiating axions [24-29]. Subsequently, domain walls are 
formed when the Hubble parameter becomes comparable to the mass of the axion H ~ m. 
Shortly after that, the condition 

A'string 

Cwall = 7 l^-Oj 

is satisfied at the time t w , where o" wa n = 9.23mF^ is the surface mass density of domain walls 
[see Eq. (E.7)], n s trmg(t) — 2vrF 2 ln(t/5 s ) is the mass energy of the strings per unit length, 
and 8 S ~ 1/V~\r] is the width of strings. After the time t w , the dynamics is dominated by 
the tension of domain walls. 

The fate of the string-wall networks is different between the case with A^dw = 1 an d the 
case with Afow > 1. If A^w = 1> domain walls bounded by strings are unstable, and they 
disappear immediately after the formation [30]. On the other hand, if A^dw > 1, the string- 
wall networks are stable, and they eventually dominate the energy density of the universe. If 
it occurs at sufficiently late time in the universe, it may conflict with standard cosmology [12]. 

2.2 Domain wall problem and its solution 

Now, let us look closer at the model with Afow > 1- After the time t w given by Eq. (2.5), 
domain walls are straightened by their tension force up to the horizon scale. In the past 
numerical studies [31], it was found that these networks of domain walls bounded by strings 
evolve into the scaling regime, which is similar to the property of the networks of strings. 
In this regime the typical length scales of defects, such as the wall curvature radius and the 
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distance of two neighboring walls, are given by the Hubble radius H 1 = 2t. Then, the 
energy density of domain walls behaves as 

Avail - ( 2 -6) 

Since the dilution of the energy density p wa ii is slower than that of matters ~ 1/R 3 
and radiations ~ 1/-R 4 , domain walls eventually dominate the energy density of the universe. 
This occurs at the time 

lb7rGcj wa ii 

where G is the Newton's gravitational constant. In other words, the wall domination occurs 
when the temperature of the universe becomes less than T ~ (D(lO)keV [20]. 

In order to avoid this overclosure problem, Sikivie [9] phenomenologically introduced a 

term 

«5V = -S»7 3 (*e- w + h.c.) (2.8) 

in the potential (2.2), where H is a dimensionless parameter which is assumed to be much 
less than unity 1 . This term, called a bias [13], explicitly breaks the discrete symmetry and 
lifts degenerate vacua represented by Eq. (2.3). If this kind of term exists in the potential, it 
affects as a pressure on walls, which eventually annihilate them. We can naively estimate the 
lifetime of the networks tdec by the ratio between the surface mass density of domain walls 
<7 wa n and the difference of the potential energy between two neighboring vacua [20] , 

tdcc ^I;! 1 ; • (2.9) 

Requiring that the decay of walls occurs before the wall domination (i.e. idee < ^wd), we 
obtain a lower bound on H 

-60 w at-3 ' m 



= >7.2xlO-»xiV- [j-^j . (2.10) 

In our previous study [20], we found that string- wall networks actually disappears in 
the time scale given by Eq. (2.9) and pointed out that significant amounts of gravitational 
waves might be produced by these defect networks. We also found that the cosmological 
abundance of cold axions produced by string-wall networks gives another constraint on the 
model parameters. In this work, we are going to give a concrete analysis of gravitational 
waves and axion radiations from these string-wall networks by numerically computing the 
spectrum of these radiations. 



3 Numerical simulations 



We performed field-theoretic lattice simulations with the model of PQ scalar field <3? whose 
potential given by Eq. (A.l). We use the three dimensional homogeneous grid with the size 
iV 3 = 512 3 . In the numerical studies, we normalize the dimensionful quantities in the unit 
of rj. For example, x — > xrj, <& — > $>/rj, etc. We choose the comoving size of the simulation 
box as L = 80 in this unit. Then, the spatial resolution is estimated as A = L/N ~ 0.16. 

1 The parameter 3 corresponds to what we called as £ in our previous study [20]. Here we do not use this 
old notation in order to avoid the confusion with the scaling parameter defined by Eq. (3.4). 
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We give the time evolution in terms of the conformal time r rather than the cosmic time t. 
We set the initial time as n = 2.0 and the final time as 77 = 40.0. The dynamical range 
of the simulation is therefore tj/ti = 20. We normalize the scale factor at the initial time 
R(ri) = 1. The other parameters are summarized in table 1. The technical details on the 
analysis method are described in Appendices A, B, C, and D. 



Table 1. Parameters used in numerical simulations. 



Grid size (N) 512 

Box size (L) 80 
Total number of time steps 1900 

Time interval (At) 0.02 

A 0.1 

m 0.1 

iV DW varying 

Initial time (rj) 2.0 

Final time (r/) 40.0 



With the parameters shown in table 1, the physical scale of the Hubble radius H , the 
width of the wall 5 W = m _1 , and the width of the string 5 S = A" 1 / 2 divided by the lattice 
spacing Ax p h ys = R(t)Ax = R(t)L/N are given by 

Axphys L Axphys Lm \t / Axphys LA 1 /"' V r / 

At the end of the simulation 77 = 40, these ratios become -ff _1 /Ax p hys — 256 < N, 
(5^/Axphys — 3.2, and (5 S /Ax p h ys — 1.01. Therefore, these length scales are marginally 
resolvable at the final time of the simulation. 

We naively guess that domain walls formed in the model with A^dw = 2 have a common 
property with that formed in the model of real scalar field with spontaneous breaking of a 
discrete Z2 symmetry. Based on this observation, we also consider the real scalar field theory 
with the Lagrangian density given by 

£ = \d r 4»'4> - V(4>), (3.2) 

-§*")'. M 

where <j> is the real scalar field. In Appendix E, we show that this model mimics the axionic 
model with iVow = 2, in the sense that the surface mass density and the width of domain 
walls are same. 

It should be emphasized that we do not perform the simulation including the bias 
term (2.8). This means that we do not simulate the annihilation process of topological de- 
fects. In order to simulate the annihilation process, much larger dynamical range is required, 
since each of Adw vacua is not homogeneously percolated in the simulation with the small 
dynamical range [20]. In our previous study [20], we could improve dynamical range by 
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performing two-dimensional lattice simulations where the size of data is given by iV 2 . In the 
present work, our aim is to investigate the spectrum of axions and gravitational waves pro- 
duced by defect networks, while it is impossible to calculate the gravitational waves properly 
in the two-dimensional field theoretic simulations. Therefore, we just calculate the spectra in 
the onset of the scaling regime, and extrapolate the results for the case in which string-wall 
networks survive for a long time. Then we use the results of two dimensional simulations of 
the annihilation process of the networks [20] when we estimate the decay time of them [see 
Eq. (4.16)]. 

3.1 Evolution of string-wall networks 

Now let us show the results of numerical simulations. We show the visualization of one real- 
ization of the simulation in Fig. 1. Using the identification scheme described in Appendix B, 
we confirmed that iVrjw domain walls attached to strings are formed at late time of the 
simulations. The structure of defect networks becomes more complicated when we increase 
the number A^dw- Figure 2 shows the spatial distribution of the phase of the scalar field 9 
in the model with Afow = 3 at selected time slices. Since we give the initial configuration of 
the scalar field as Gaussian random amplitude (see Appendix A. 3), the scalar field randomly 
oscillates around the minima |<1>| 2 = rf of the potential, given by the first term of Eq. (2.2). 
Strings are formed when this initial random configuration becomes relaxed. Subsequently, 
the value of the phase of the scalar field 6 is separated into Afow = 3 domains. Domain walls 
are located around the boundaries of these domains. 

We note that the initial stage of simulations does not correctly follow the history of the 
universe, since domain walls are formed immediately after the formation of strings. As we 
noted in Sec. 2.1, in the realistic case, strings are formed much earlier than the formation of 
domain walls. This is due to the hierarchy between two different scales, Aqcd/^ ~ 10 -11 . 
On the other hand, in our choice of the parameters in numerical simulations, this ratio 
becomes as small as Aqcd/-F<j ~ ~ O(0. 1) (recall that we used the unit rj = 1). Because 
of the limitation of the dynamical range of the numerical simulations, we cannot choose this 
ratio as an arbitrary small value. Therefore, we concentrate on the evolution at late times 
where the string-wall networks evolve along to the scaling regime, and ignore the effect of 
the initial relaxation regime. We will discuss this point again in Sees. 3.2 and 3.3. 

We define the length parameter £ of strings 

£ = ^*™£ t 2 > (3 A) 

/^string 

and the area parameter A of domain walls 

A = ^t. (3.5) 

Cwall 

These parameters take a constant value when the defect networks enter the scaling regime 
(i.e. /String ~ 1/i 2 and p wa n ~ 1/t). Henceforth we refer to these parameters as the "scaling 
parameters". We plot the time evolution of these parameters in Fig. 3. From this figure, 
we see that the defect networks are in scaling regime in the time scale r > 15. The length 
parameter £ does not strongly depend on A r DWi and its value is consistent with £ ~ 0.5- 
1.3 obtained in previous studies [11, 28, 29, 32] On the other hand, the value of the area 
parameter A increases for large Afow- This fact agrees with the intuitive argument that 
the number of domain walls increases proportionally to Adw if the number of strings per 
simulation box is unchanged. 
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(a) AW = 2 



(b) AW = 3 




(e) AW = 6 



Figure 1. Visualization of the simulations with (a) A^dw = 2, (b) AWv = 3, (c) AWv = 4, (d) 
^Vdw = 5, and (c) A^dw = 6. In this figure, we take the box size as L = 60 and A" = 256, which 
is smaller than that shown in table 1. Each figure shows the spatial configurations of topological 
defects at the time r = 25. The white lines correspond to the position of the core of strings, which 
is identified by using the method described in Appendix B.l. A^dw domain walls are represented by 
surfaces with various colors, which are identified by using the method described in Appendix B.2. 

The right panel of Fig. 3 also implies that the area of domain walls in the real scalar field 
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(c) r = 19 (d) r = 25 



Figure 2. The distribution of the phase of the scalar field 9 on the two-dimensional slice of the 
simulation box. Each figure shows the visualization at different times: (a) r = 7, (b) r = 13, (c) 
r = 19, and (d) r = 25. In these figures, we used the same data that are used to visualize the result 
with TVdw = 3 in figure 1. The value of 9 varies from — tt (blue) to it (red). At late times, the value 
of 9 is separated into three domains represented by blue, red, and green region. Domain walls are 
located around the boundary of these three regions, 9 = 7r/3, 9 = — 7r/3, and 9 = ±ir. Strings, which 
are represented by white lines, pass through the point where three regions meet each other. 

theory given by Eqs. (3.2) and (3.3) evolves similarly to the axionic model with iVow = 2. 
This means that the behavior of the networks for the case with iVow = 2 is similar to that 
produced by the model with Z<i symmetry, except that they contain strings on the surface of 
walls. We note that the bump-like behavior of A in the real scalar field model at the initial 
stage of the simulation should be regarded as a numerical fake. We identify the position 
of domain walls as the point on which scalar field changes its sign in the model with the 
real scalar field, which is different from the identification scheme used for the model with 
the complex scalar field (see Appendix B.2). These different identification schemes cause the 
different behavior of A at the initial stage, where domain walls do not relax into the scaling 
regime. 
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Figure 3. Time evolution of the length parameter £ (left panel) and the area parameter A (right 
panel) for various values of -ZVdw- The points denoted as "real" in the right panel show the evolution 
of A in the model with real scalar field defined by Eqs. (3.2) and (3.3). 

3.2 Production of axions 

Prom data of the scalar field obtained by numerical simulations, we compute the power 
spectrum P(k) of axions radiated by string-wall networks, which is defined by 

i<d(t,k)*d(t,k')> = ^V)(k-k')P(M), (3-6) 

where d(i, k) is the Fourier component of the time derivative of the axion field, and (. . . ) 
represents an ensemble average. Here, we mask the position where topological defects exist, 
and extract the pure radiation component in the axion field produced by the defect networks. 
See Appendix C for details on the analysis method. 

Figure 4 shows the results for P(k,tf) for various values of A^dw evaluated at the final 
time of the simulation. We find that the spectrum has a peak around the momentum scale 
determined by the mass of the axion, k/R(tf) ~ m, or k ~ R(tf)m ~ 2. This implies that 
most of the axions are produced by the self-interaction of domain walls whose width is given 
by the inverse of the mass of the axion ~ m . However, this peak does not seem to fall off, 
but makes a tail at higher momenta. This tail-like feature was also found in another numerical 
study [11] in which we investigated the spectrum of axions produced by the annihilation of 
-^dw = 1 domain walls bounded by strings. The appearance of this high-frequency tail can 
be understood in terms of the conjecture made in Refs. [26, 27, 33, 34]. The strings whose 
radius smaller than the horizon size quickly decay due to the tension of domain walls in the 
time scale comparable to (or less than) the Hubble time. This fast process does not change 
much the energy spectrum of the axion field composing the string configuration, which has 
a form dE/dk ~ 1/k with a high-momentum cutoff of order d^ 1 (the width of strings) and a 
low-momentum cutoff of order t~ l (the Hubble radius) [26]. This 1/k behavior makes a tail 
in the spectrum of radiated axions, as we see in Fig. 4. This argument cannot be applied for 
the spectrum of axions produced by global strings [28, 29], which has a sharp peak at the 
horizon scale since in the absence of domain walls the strings lose their energy in the time 
scale much longer than the Hubble time. The high-momentum cutoff which is given by the 
width of the string ~ J" 1 is not clearly shown in Fig. 4 because of the poor resolution at the 
final time of the simulation. 
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Figure 4. The spectra of axions produced by string- wall networks for various values of Njyw- These 
spectra are obtained by using the masking analysis method described in Appendix C. The spectra 
shown in this figure are the difference of the power spectrum evaluated at two different time steps T\ 
and T2 [see eq. (C.12)]. Here, we choose T\ = 14 and T2 = 40. 



Note that the amplitude of the high-frequency tail becomes small if we increase the 
domain wall number JVdw- I n the case of large JVqwj the decay of loops of strings with small 
curvature radius might be "frustrated", because of the tension of domain walls acts on the 
strings from various directions. Then the radiation of hard axions with spectrum ~ 1/k is 
suppressed, compared with the case of small -/Vtjw- 

The above observation can be confirmed by estimating the energy of radiated axions. 
Roughly speaking, the product P{k)k represents the fraction of the energy density of axions 
with comoving wavenumber k against the total energy of axions. Hence the ratio of the 
energy between hard axions which contribute to the tail of the power spectrum and the 
energy of axions which contribute to the peak of the power spectrum can be estimated 
as £WW ~ [^(fchard)*Wd]/[P(W)W] ~ P(10~ 6 ) x O(10)]/[O(10~ 4 ) x 0(1)] ~ 
10 _1 . On the other hand, the ratio between the energy of domain walls and the energy of 
strings at the final time of the simulations becomes £ s tring/£waii ~ [i/A i string]/[^/0'wail] ~ 10 . 
Therefore, strings have 10% of the energy of domain walls at the final time, which contribute 
to the hard component ~ 1/k of the power spectrum of radiated axions. Since the fraction 
-^string /-E'waii decreases with time, we expect that the amplitude of the high-frequency tail 
would be negligible compared with the height of the peak of the spectrum if we follow the 
evolution for much larger dynamical ranges. 

In the spectrum shown in Fig. 4, we subtract the contribution of axions produced 
before a pivot time t\ [see Eq. (C.12)], since the axions produced at initial relaxation regime 
contaminate the spectrum of axions radiated from the string-wall networks. Here, t\ should 
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Figure 5. The power spectrum of axions evaluated at the time = 40 for the case TVdw = 3. The 
red plot shows the spectrum evaluated with masking the contamination of the core of topological 
defects and subtracting the contribution of axions produced at initial stage (r < 14). The green plot 
shows the spectrum obtained by the same model without masking the contamination of defects. The 
blue plot shows the spectrum evaluated without masking and subtraction. 



be chosen as the time of the onset of the scaling regime. In the numerical study, we use 
7~i = 14 as the subtraction time, at which the system begins to follow the scaling regime. 
Figure 5 shows the effect of this subtraction procedure. We see that the position of the peak 
is falsified, and the amplitude of the spectrum is overestimated if we do not perform the 
subtraction. Figure 5 also shows that the spectrum is overestimated at higher momentum 
scale when we do not mask the data near the position of topological defects. 

Since the population of the axions radiated by string-wall networks is dominated by 
low- momentum modes which has the frequency comparable to the mass of the axion m, we 
expect that the mean energy of radiated axions is determined by the parameter m. Here, 
using the result of P(k,t2), we compute the mean comoving momentum of radiated axions 



fc(i 2 



and its ratio to the axion mass 



k(t 2 )/R(h 



(3.7) 



(3i 



The results are shown in table 2. We see that the value of the mean momentum is not 
strongly depend on the number A?dw- However, the uncertainty becomes large for the case 
with large -/Vdw- This is because a large portion of the simulation box is masked when we 
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evaluate the power spectrum, which causes a large systematic errors when we evaluate the 
spectrum in the large scale 2 , as shown in Fig. 4. 

Table 2. The values of mean comoving momentum k and parameter e a defined in Eqs. (3.7) and 
(3.8) for various values of -/Vdw- 



AW 


k 


e a 


2 


2.38±0.28 


1.19±0.14 


3 


2.38±0.25 


1.19±0.12 


4 


2.63±0.29 


1.31±0.15 


5 


2.97±0.45 


1.49±0.22 


6 


3.07±0.42 


1.54±0.21 



The claim is that the parameter defined by Eq. (3.8) does not depend on the time <2. 
To confirm it, we calculate the parameters k and e a at different time steps, while fixing the 
reference time t\ = 14 in Eq. (C.12). The results are shown in table 3. We find that the value 
of k changes proportionally to Rfo)- On the other hand, the ratio e a between the physical 
momentum k/R and the axion mass m remains constant. Therefore, we use the parameter 
e a as a constant of 0(1) when we perform the analytic investigations in Sec. 4. 

Table 3. The values of k and e a evaluated at different times for the case with JVdw = 3. 



T2 


k 


e a 


24 


1.42±0.22 


1.19±0.18 


28 


1.62±0.13 


1.15±0.09 


32 


1.88±0.19 


1.18±0.12 


36 


2.04±0.22 


1.13±0.12 


40 


2.38±0.25 


1.19±0.12 



We close this subsection by showing Fig. 6 the time evolution of the total energy density 
of axions and topological defects in the simulation box. From this figure, we see that the 
domain wall networks evolve as p W aii ~ 1/t, or p wa ii ~ 1/V 2 in conformal time. This behavior 
is expected from the property of the scaling solution, given by Eq. (2.6). We also find that 
the energy density of axions radiated by string-wall networks behaves similarly to that of 
topological defects at the late time. We will confirm this behavior in Sec. 4.1. 

3.3 Production of gravitational waves 

Next, we investigate the production of gravitational waves from string- wall networks. The 
spectrum of gravitational waves is represented by the quantity 



1 dp gw (k,t ) 
p c (t) dink 



"gw(fc,*) = T^ ( 3 - 9 ) 



We find that there is a systematic error in our masking analysis method. See Appendix C for details. 
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Figure 6. Time evolution of energy densities for various values of -ZVdw- Pdefects represents the 
energy density of the string- wall networks which is evaluated by summing up the data of scalar field 
placed near the defects whose position is identified by using the algorithm described in Appendix B. 
Pa, kin, Pa.grad, and p a ,pot are the kinetic, gradient, and potential energy density of radiated axions, 
respectively. p a ,tot is the sum of p a ,kin, Pa,grad, and p a , po t- These energy densities of axions are 
evaluated by subtracting the data of scalar field which contribute to Pdefects- 

where p c {t) is the critical density of the universe at the time t, and p gw is the energy density 
of gravitational waves. We compute figw using the data of the scalar field $ obtained from 
the numerical simulation. The method of the calculation is summarized in Appendix D. 

Figure 7 shows the results for S7 gw (/c, t) for various values of A'dw- The slope of the spec- 
trum changes at two characteristic scales. One is the Hubble radius {k ~ 2irR(Tf)H(Tf) ~ 
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0.1) and another is the width of the domain wall (k ~ R{Tf)m ~ 2). This result agrees 
with the general argument made by two of the present authors (M. K. and K. S.) in [35]. 
However, the amplitude of S7 gw at high momenta is enhanced for the case with large A^ow- 
This can be interpreted as follows. Since the number of domain walls is large in the model 
with large Adw, domain walls are "frustrated", and there are many small scale configura- 
tions whose correlation length is shorter than the Hubble radius. Therefore, the small scale 
anisotropy becomes large in the model with large domain wall number A?dw- We also plot 
the spectrum of gravitational waves produced by domain walls in the real scalar field theory 
given by Eqs. (3.2) and (3.3). Although the evolution of domain wall networks shown in 
Fig. 3 is similar between the model with Afow = 2 and the model with real scalar field, the 
form of the spectrum of gravitational waves is slightly different between these two models at 
high momenta. This might be caused by the different form of the effective potential which 
determines the small scale structure of the defects. 
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Figure 7. The spectra of gravitational waves for various values of Nuw- The plot denoted as "real" 
shows the spectrum of gravitational waves produced in the model with real scalar field defined by 
Eqs. (3.2) and (3.3). In this figure, we normalized the vertical axis by the amplitude at the peak 
momentum. 

As we discussed in previous subsections, at the early times (r < 15) topological defects 
do not relax into the scaling regime. However, the nonlinear dynamics of scalar fields during 
this initial stage produces a "burst" of gravitational waves, which may contaminate the 
spectrum of gravitational waves produced at late times where the defects networks evolve 
along to the scaling solution. We subtract this contribution with the same manner that we 
used for the subtraction of axions produced at the initial stage [see Eq. (D.19)]. Figure 8 
shows the time evolution of the spectrum of gravitational waves, where we compare the effect 
of this subtraction procedure. We see that the form of the spectrum is different at early times 
if we make a subtraction. In particular, there is a bump at the scale given by the Hubble 
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radius k/R ~ H, and a plateau which falls off at the scale given by k/R ~ m. This form is 
different from what we found in our previous studies [35, 36], where a nearly "flat" spectrum 
extends between two characteristic scales. It seems that this difference in the form of spectra 
is caused by the contamination of radiations produced at the initial relaxation regime. We 
also note that the dynamical range of the past numerical studies (tj/tj ~ 12) [35, 36] is 
shorter than that of the present studies (t//tj = 20). Therefore, the bump-like feature is 
less apparent in the past numerical studies with shorter dynamical ranges. The form of the 
spectrum at the final time of the simulation is similar between the result with subtraction 
and that without subtraction, which indicates that contaminations of the initial stage are 
diluted due to the cosmic expansion. 



N DW = 3 without subtraction N DW - 3 with subtraction 




comovina wavenumber k comovina wavenumber k 



Figure 8. The time evolution of the spectrum of gravitational waves for the case with ./Vow = 3. 
Here, we plot the function Sk(r) defined by Eq. (D.17) (left panel) and AS'fc(r) defined by Eq. (D.19) 
(right panel), which is proportional to f2 gw in the radiation-dominated background. In the right panel, 
we subtract the contribution of gravitational waves produced before the time r g = 15. The spectra 
are shown from the conformal time r = 20 (pink) to r = 40 (green) with the interval At = 5. 



The spectrum shown in Fig. 7 is normalized in terms of the peak amplitude in order 
to give the relative form of the spectrum between different theoretical parameter JVdw. On 
the other hand, the total amplitude of gravitational waves is determined by the following 
arguments. Since the spectrum of gravitational waves has a peak at the momentum which is 
given by the Hubble scale, we regard that the most of the gravitational waves are generated 
at the length scale ~ t. Suppose that the radiation of gravitational waves occurs with the 
time scale comparable to the Hubble time. The power of gravitational waves radiated from 
domain walls is estimated by using the quadrupole formula [37] 

P ~ GQijQij, (3.10) 

where Qij ~ Mdw^ 2 is the quadrupole moment, and Mow ~ Pwaii^ 3 ~ ^waii-^ 2 is the 
mass energy of domain walls, where we used Eq. (3.5). The energy of gravitational waves is 
estimated as 

E gw ~ Pt ~ GA 2 a 2 wall t 3 . (3.11) 
We see that the energy density of gravitational waves does not depend on the time t 

Pgw ~ £ g wA 3 ~ GA 2 a 2 wall . (3.12) 
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Now, we define the efficiency e gw of gravitational waves 



(sim) 



(3.13) 



where Pgw™^ is the total energy density of gravitational waves computed from numerical 
simulations. We plot the time evolution of the parameter e gw in Fig. 9. Although the value 
of e gw does not remain in constant exactly, it seems to converge into a universal value of 
0(5-10). 
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Figure 9. The time evolution of the efficiency parameter e gw defined by Eq. (3.13) for various values 
of -/Vdw- The points denoted as "real" in the right panel show the evolution of e gw in the model 
with real scalar field defined by Eqs. (3.2) and (3.3). In this figure, we subtract the contribution of 
gravitational waves produced before the time r g = 15. 



We note that the time dependence of e gw is strongly affected by the choice of the time T g 
at which we subtract the contribution of gravitational waves produced in the initial relaxation 
regime [see Eq. (D.19)]. This is shown in Fig. 10, where we plot the time evolution of e gw for 

various values of r g . If we choose r g at an early time of the simulation, contains the 

contribution of gravitational waves produced during the initial relaxation regime. Hence the 
behavior of /Ogw"^ might be different from what we expect from Eq. (3.13), where we assumed 
that the energy density of domain walls is given by the scaling formula (3.5). This additional 
contribution is not negligible even at the late time of the simulations because the dynamical 
range of the numerical simulation is short. This constancy of e gw should be tested in future 
numerical studies with improved dynamical ranges. Since the time variation of e gw is small 
for large value of T g , as shown in the plot of r g = 25 in Fig. 10, we assume that e gw takes 
a constant value when string-wall networks enter into the scaling regime, and use the value 
£gw — 5 when we estimate the abundance of gravitational waves in the next section. 
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Figure 10. The time evolution of e gw and its dependence on the choice of T g . In this figure, we plot 
the result of the model with TVdw = 3. 



4 Observational constraints 

In this section, we give constraints on the parameters of axion models with long-lived domain 
walls, using the results of numerical studies described in the previous section. First, we 
calculate the relic abundance of cold axions and gravitational waves produced by string-wall 
networks in Sec. 4.1. Then we briefly review some experimental constraints in Sees. 4.2 
and 4.3. We combine these constraints, and discuss their implications on the models with 
Ndw > 1 in Sec. 4.4. We find that the overabundance of the cold axions produced by domain 
walls gives a stringent bound on the model parameters, which excludes the most part of 
models, except for tiny loopholes. Finally we comment on some exceptions in Sec. 4.5. 

4.1 Axion cold dark matter abundance 

Here we assume that the energy of strings becomes negligible compared with that of domain 
walls, and axions and gravitational waves are only produced by domain walls. In this approx- 
imation, the evolution of energy densities for domain walls, axions, and gravitational waves 
is described by the following coupled equations 



(4.1) 



emission 



dt r dt 

dp a ou i dp w ^ a 

-dT = - 3Hpa + ^r> (42) 

dfgW _ ATT n | dp w ^g 

__-4iJ Pgw + ^^, (4.3) 
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where 



dt 



emission dt dt 



is the rate of radiation from domain walls. Following from Eqs. (3.5) and (3.13), we put an 
ansatz on the behavior of the energy densities of domain walls and gravitational waves 

Pwaii = A^j^-, and p gw = e gw GA 2 a^ an , (4.5) 

where A and e gw are constants. Substituting these expressions to the above equations, we 
obtain 



dp. 



wall 



dt 



= A^, (4.6) 

emission 

dt ~ gw t 



2e gw —■ (4.7) 



From Eqs. (4.4), (4.6), and (4.7), we find 



dp w ^a _ . *7 wa n _ ^^ 2 °" wall 

dt ~ A 2t* gw t 



Note that the sign of the right hand side of this equation becomes negative for t > (Gcr wa ii) . 
This is not problematic, since domain walls have to collapse before the time t ~ (Gcr wa u) _1 
otherwise they overdose the universe [see Eq. (2.7)]. 

Let us define the energy and number of axions per comoving box 

E a (t) = R\t)p a (t), (4.9) 
N a (t) = R 3 (t)n a (t), (4.10) 

where n a {t) is the number density of axions at the time t. Integrating Eqs. (4.2) and (4.8), 
we find 



E a {t)= f dt'a 3 (t') 



AOwatt ,_2 r, n jl A2+-1 

— - — t - 2e gw Go wall .4. t 



(4.11) 



where t r is the time when domain walls begin to radiate axions, which might be chosen as 
the time of the QCD phase transition. In the regime t > t r , it reduces to 

E a {t) ~ R\t) (A^ - le gw Gal all A^ . (4.12) 

For sufficiently early times t <C (Go^an) -1 , the first term of Eq. (4.12) dominates over the 
second term, which leads to the behavior p a oc 1/t. This agrees with what we observed in 
numerical simulations, shown in Fig. 6. The ratio between E a {t) and N a (t) can be determined 
by the result of numerical simulations 

f§ - < 4 - 13 > 

where e a is defined by Eq. (3.8). Suppose that domain walls disappear at the time idee and 
the radiation of axions and gravitational waves is terminated at that time. In this case, the 
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energy densities of axions and gravitational waves at the present time to are given by 



Pa,dec(*o) 



Pgw(*o) 



mN a (t dcc ) 
R 3 (t ) 

\ R(to) 

V R{to) 



A 



Q"wall 

^dcc 



-c Crr 2 A 2 
t gw LxCT wall >V 



egwGA 2 ^ 



wall' 



(4.14) 
(4.15) 



where the subscript "dec" in the left hand side of Eq. (4.14) represents that the axions 
produced by the decay of domain walls. 

The annihilation of domain walls occurs in the time scale estimated by Eq. (2.9) 



*dc 



■/?? 



a 



N DW Ei] 2 ' 



(4.16) 



The numerical coefficient can be fixed to a ~ 18, according to our previous study [20]. Note 
that the relation 



R(t 



dec) 



R(to) 



5.45 x 10~ 9 x N, 



-3/2 
DW 



10 10 GeV\ 3/2 /10~ 58 ^ 1/2 



where we used the expression for the mass of the axion 

m = 6 x 10 _4 eV 



10 10 GeV \ 

Fa )' 



(4.17) 



(4.18) 



Substituting Eqs. (4.16) and (4.17) into Eqs. (4.14) and (4.15), we obtain the density param- 
eter of axions and gravitational waves 



^a,dec(*o)^ 2 



Pa,dcc{to) 
Pc(t )/h 2 

1.01 x 10 3 x 



A 



-3/2 fW 



DW 



58 \ 1/2 



1 - 5.35 x 10~ 3 x e gw AN^ 



10 



-58 



10 10 GeV \ 1/2 
Fa J 



10 10 GeV 

F a 



(4.19) 



^gw(^o)^- 2 



Pgwjtp) 
Pc(t )/h 2 



2.24 x 10~ 8 x e gw A 2 N D ^ 



10 



-58 \ 2 



HTOeVY 

Fa ) 



(4.20) 



where h is the reduced Hubble parameter today: Hq = 100ft,km-sec -1 Mpc -1 . We can fix the 
numerical values of A, e a and e gw from the result of numerical simulations. Then, the relic 
abundance of axions and gravitational waves is determined by three theoretical paramters: 
A^dWj S, and F a . The second term in the square bracket in Eq. (4.19) represents the effect of 
gravitational radiation, which is negligible at early times but becomes relevant when domain 
walls survive for a long time. In particular, when domain wall survives enough time such 
that they overdose the universe, tdec ~ (Gcr wa n) _1 , the contribution of the second term in 



- 20 - 



the square bracket in Eq. (4.19) becomes 0(1). In such a case, the gravitational field around 
walls cannot be regarded as a small perturbation [38], and we must use the full general 
relativistic treatment, which is out of the scope of this paper. Hence, the expression (4.19) 
is valid only for the time t <C twD ~ (Ga^n)' 1 . This condition is automatically satisfied as 
long as we consider the parameter region where domain walls do not overdose the universe. 

We note that there are other contributions to the axion cold dark matter abundance, 
in addition to fi a ,dec- O ne is the coherent oscillation of the homogeneous axion field (zero 
mode) [7] and another is the production of axions from global strings [24] . Their contributions 
are estimated as [11] 

^a,o(io)/i 2 - 0.6 x 

^a,string(io)/i 2 - 4.0 X 

where the subscripts "0" and "string" represent that they are the contributions of zero modes 
and global strings, respectively. In the above expressions, we fixed the QCD scale such that 
Aqcd = 400MeV. The total abundance of cold dark matter axions is given by the sum of 
Eqs. (4.19), (4.21) and (4.22), 

— £l a ()h + ^a, string 

h 2 + n a:dcc h 2 . (4.23) 

We require that £l a h 2 should not exceed the observed value of the abundance of cold dark 
matter Ocdm^ 2 = 0.11 [39]. This gives an upper bound on F a and a lower bound on 5. 

4.2 Neutron electric dipole moment 

Next we consider the effect of the E term on the degree of strong CP-violation, following 
Refs. [17, 18]. The QCD Lagrangian contains a term 

C = 32^^^' (424) 

where Q a ^ v is the gluon field strength, is the dual of it, and is a dimensionless pa- 
rameter. This term contributes to the neutron electric dipole moment (NEDM) with the 
amplitude d n = 4.5 x 10~ 15 #ecm [3]. The recent experimental bound on the NEDM gives 
\d n \ < 2.9 x 10~ 26 ecm [40], which requires 

6 <0.7 x KT 11 . (4.25) 

The original idea of Peccei and Quinn is to set 8 = by introducing a symmetry [2]. How- 
ever, if the additional term (2.8) exists in the effective potential, it shifts the CP-conserving 
minimum with a magnitude controlled by the parameters H and 5. 

Substituting the parametrization $ = r]e ia l r] into the effective potential given by Eqs. (2.2) 
and (2.8), and expanding for small a/77, we obtain the effective potential for a 

V(a) a \m 2 phys a 2 - m 2 phys F a 6a, (4.26) 

where m p h ys is the effective mass of the axion 

m 2 phys = m 2 + ml = m 2 + 2EN^ W F 2 cos 6, (4.27) 



F a 



10 12 GeV 

F a 

10 12 GeV 



1.19 



1.19 



(4.21) 
(4.22) 
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and 8 is the shifted minimum 







(a) 

F a 



2HiY3 w F a 2 sm«5 



(4.28) 



m 2 + 2~AT 2 w F 2 cos<5' 



Requiring that this should not exceed the experimental bound (4.25), we obtain an upper 
bound for H. 

4.3 Astrophysical bounds 

It is known that various astrophysical observations give a lower bound on F a . Here we just 
summarize the constraints obtained so far. See [41, 42] for more comprehensive reviews. 

First of all, the emission of axions in the sun gives a constraint on the axion-photon 
coupling, g arr = (a em /27ri ? a )c a77 , where a em is the fine-structure constant, and c a77 is a 
numerical factor of 0(1). The helioseismological consideration on the anomalous solar en- 
ergy loss due to the axion emission gives a bound <7 a77 < 10 _9 GeV _1 [43]. There are also 
direct experimental searches. The Tokyo Axion Helioscope [44] gives a bound g a77 < (5.6- 
13.4) x Kr 10 GeV _1 for 0.84eV < m < l.OOeV. Phase I of the CERN Axion Solar Telescope 
(CAST) [45] gives g ail < 8.8 x 10" 11 GeV" 1 for m < 0.02eV. It is improved in Phase II [46] 
as da-y-y < 2.2 x 10 _10 GeV -1 for m < 0.4eV. The sensitivity is expected to be improved up 
to <? a77 < few x 10~ 12 GeV _1 in the next generation helioscopes [47]. 

The observation of globular clusters gives another bound. The helium-burning life- 
times of horizontal branch (HB) stars give a bound for axion-photon coupling g a ^ n < 0.6 x 
10 _10 GeV~ 1 [48]. Furthermore, the delay of helium ignition in red-giant branch (RGB) stars 
due to the axion cooling gives g aee < 2.5 x 10~ 13 [49], where g aee = C e m e /F a is the axion- 
electron coupling, m e is the mass of the electron, C e = cos 2 /J/3, and tan f3 is the ratio of two 
Higgs vacuum expectation values. 

The axion-electron coupling is also constrained by the observation of white-dwarfs. The 
cooling time of white-dwarfs due to the axion emission gives a bound g aee < 4 x 10" 13 [50]. 
Recently, it is reported that the fitting of the luminosity function of white-dwarfs and the 
pulsation period of a ZZ Ceti star is improved due to the axion cooling, which implies the 
axion-electron coupling g aee ~(0.6-1.7) x 10~ 13 [51] or g aee ~(0.8-2.8) x 10~ 13 [52]. 

Finally, the axion-nucleon coupling is constrained from the energy loss rate of the super- 
nova 1987A [41]. This gives the strongest astrophysical bound on the axion decay constant 
compared with other bounds described above, 



Hence in the following we take it as a lower bound on F a when we combine it with another 
observational constraints. Note that the estimation of the axion emission rate suffers from 
various numerical uncertainties which may modify the bound by a factor of 0(1) [41]. 

4.4 Implication for models with Nbw > 1 

In Fig. 11, we plot the observational constraints described in the previous subsections. Based 
on the estimation given by Eq. (4.23), we find that the overabundance of cold axions gives 
a stringent lower bound on the bias parameter, 5 > 10~ 50 -10~ 52 . This lower bound is much 
stronger than that obtained from the overclosure of domain walls (2.10). Combined with the 
NEDM constraint, it completely excludes the parameter region if we assume that the phase 
5 of the bias term (2.8) is 0(1). This bound might be weakened if we assume the extremely 



F a > 4 x 10 



> 8 GeV. 



(4.29) 



- 22 - 



10' 



r45 




10' 



Q a > ^CDM 




10' 



,8 



10* 10 



10 



,11 



10 



,12 



F a [GeV] 



Figure 11. The various observational constraints in the parameter space of F a and 3. The green 
dashed-line represents the parameter region estimated by Eq. (4.23) where f2 a = Ocdm is satisfied, 
and the region below this line is excluded since the relic abundance of axions exceeds the cold dark 
matter abundance observed today. The vertical dottcd-linc represents the bound given by Eq. (4.29) 
which comes from the observation of supernova 1987A. The red solid-lines represent the NEDM bound 
given by Eqs. (4.25) and (4.28) for 8=1, 10 -4 , and 10 -8 . The region above these lines is excluded 
since it leads to an experimentally unacceptable amount of CP-violation. The pink line represents 
the NEDM bound for the critical value 8 = 1.1 x 10~ 2 . There are still allowed regions if the value of 
8 is smaller than this critical value. The blue dotted-lines represent the parameter regions on which 
the amplitude of gravitational waves, given by Eq. (4.20), becomes f2 gw /i 2 = 1CP 20 and 10 -22 . In this 
figure, we fixed other parameters as -/Vdw = 6, e a = 1.5, A = 2.6, and e gw = 5. 

small value of 5, but this assumption spoils the genius of the original PQ solution to the 
strong CP problem. We must fine-tune 5 in order to solve the fine-tuning problem of 8. 

We find that the amplitude of gravitational waves produced by domain walls is weak 
Q gw h 2 < 10 -20 in the parameter space considered above. This is irrelevant to any observa- 
tions, since even the ultimate phase of DECIGO [53], which is a space-borne interferometer 
planned to launch in the future, will have a sensitivity Q gw h 2 ~ 1CP 18 [54]. The weakness 
of the signal of gravitational waves is due to the fact that domain walls annihilate in early 
epoch in order to avoid the overclosure of axions produced from them, and the dilution factor 
of cosmic expansion becomes large. 

The difference of this result from our previous study [20] arises due to the different 
estimation on the mean momentum of radiated axions. In the present analysis described in 
Sec. 3.2, it turns out that the radiated axions are mildly relativistic with the momentum 
comparable to their mass, k/R ~ ?72j or c a — 1 . However, in the previous study [20], we 
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assumed that the radiated axions have a mean energy uj ~ 7?n where 7 ~ 60, according to 
the result of the past numerical study 3 [34] . This leads to an underestimation of the number 
of axions produced by domain walls by a factor + 1/7 ~ 10~ 2 . In this case, the bound 
which comes from the overclosure of cold axions is weakened. The large factor of 7 originated 
in the assumption made in [34], that the radiated axions have a hard spectrum dE/dk ~ 1/k 
which leads to a logarithmic correction 7 ~ \n(^f\ri/m) ~ 60. In Sec. 3.2 we have shown 
that this contribution is negligible compared with the contribution of the peak located at the 
axion mass. 

We note that there is a parameter region around F a ~ 10 6 GeV, called the "hadronic 
axion window" [56-58], that is not excluded by observations. This occurs due to the am- 
biguity in light quark masses which leads to the cancellation in the axion-photon coupling 
g ail for KSVZ models. In such a case, the astrophysical bounds on g ail do not significantly 
constrain the value of F a . Here we do not consider this possibility, as it strongly depends 
on the quantum number assignment of the PQ-charged fermions [58]. Also, we note that in 
this case the axion becomes a candidate of hot dark matter [57] which can be constrained by 
observation of the large-scale structure [59]. 

From the results shown in Fig. 11, we conclude that the axion models with ./Vdw > 1 
are excluded if the PQ symmetry is broken after inflation. The exception is the case in which 
the value of 5 is suppressed. As we see in Fig. 11, the degree of tuning in 5 is < 1.1 x 10 -2 
in order to avoid all observational constraints. 

4.5 Scenario with extremely small 5 

Although there is no theoretical reason to consider the extremely small value of 5, it is 
possible, by accident or design of some fundamental physics, to yield 5 whose value is small 
enough to avoid the NEDM bound. Let us comment on such a situation. In this case, the 
value of S can take larger value than that we considered so far. In particular, if S becomes 
as large as 



the axion mass is not determined by QCD instanton effect, but by the correction term m* 
[see Eq. (4.27) and also Ref. [18]]. Then, the cosmological history is significantly modified. 

When the condition (4.30) is satisfied, the magnitude of the bias term (2.8) exceeds that 
of the QCD potential [the second term of Eq. (2.2)]. Hence, "domain walls" with the tension 
0"wall — are formed when the Hubble parameter becomes H ~ m*. The subsequent 

history is categorized into two possibilities, analogously to the scenarios described in Sec. 2. 
One possibility is that strings are get attached by domain walls at the time t ~ m*, and 
they quickly disappear due to the tension of walls. Let us call this scenario case I. On the 
other hand, if the structure of the term (2.8) is the form oc 3> n , where n is an integer, n 
domain walls are attached to strings, and they survive for a long time. Let us call this 
scenario case II. These string-wall networks annihilate due to the effect of a "bias" term 
which arises from the QCD instanton effect, 5V ~ Aq CD ~ m 2 F 2 . This occurs at the time 
*dcc,* ~ ^wall/A-QCD ~ 8m„/m 2 . Requiring that domain walls should disappear before the 

3 We note that there is another numerical study which indicates 7 ~ 3 [55]. Our present results support 
this estimation rather than 7 ~ 60 obtained in [34]. 




(4.30) 
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epoch of big bang nucleosynthesis (BBN), idee,* < ^bbn ~ lsec, we obtain 



2<2x10 JxW bw(jp4vJ ' (431) 

In case I, the abundance of axions is estimated in the usual way, as a sum of the 
coherent oscillation, strings, and decay of domain walls bounded by strings. The number 
of axions is fixed at t ~ i* = m~ l . At that time, the energy density of axions is given by 
Pa(t*) ~ 10 x m 2 F 2 , where the factor 10 arises from the fact that the abundance of axions 
becomes larger by an order of magnitude if we include the contribution of strings and domain 
walls [11, 29]. The present energy density becomes p a (to) ~ 10(R{t*) / R(to)) 3 rn%Fa > which 
leads 

/ \ V 4 / rp \ 5/2 

nZU D^ 2 - 1-2 x io- 5 x <(J ^ j (tp^v. 

where the subscript "nt" indicates that they are produced non-thermally. 

Axions are also produced from thermal bath of the primordial plasma. The difference 
from the usual scenario is that the axion mass is given by m* so that their relic energy density 
can be large. The number of axions is fixed when their interactions with the thermal bath 
with the process mediated by gluons freeze out. Their number at the present time is given 
by n* h (to) — 7.5cm -3 [60], which leads 

Of *> = 3.2 X 10- X AW ( I | 3? ) , (4.33) 

where the subscript "th" indicates that they are produced thermally. For sufficiently large 
value of S, the contribution of thermal component (4.33) exceeds the non-thermal conpo- 
nent (4.32). Requiring that it should not exceed the present cold dark matter abundance 
n^h 2 < 0.11, we obtain 

S < 1.2 x 10- 34 x N D l (j^^— ) 2 easel. (4.34) 

On the other hand, in case II, axions are copiously produced from long-lived domain 
walls. The relic abundance of axions produced by domain walls can be estimated by using 
the argument analogous to that leads Eq. (4.19). Then we obtain 

/ \ 3 / 4 / TP \9/ 2 

^ (caS en)^-1.2xl0-x^^j • ( 4 - 35 ) 

In this case, the thermal component does not exceeds the non-thermal component. Requiring 
that ^ case lT) h 2 < 0.11, we obtain 

H < 9.3 x 1(T 42 x AT D ^ ( 1Ql f° eV ) case II. (4.36) 

This bound is more severe than that obtained from the lifetime of domain walls, given by 
Eq. (4.31). 
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However, if S is much larger than the values considered above, axions become heavy and 
no longer stable. Their dominant decay channel is a decay into two photons. The lifetime is 
given by 

„ n 64vr 

u = r 



7 mhl-y-y 

= 1.1 x 10«sec x ( JL, ) " 3/2 (^^) , (4-37) 

where T 7 is the decay late of the process 0—7-77, an d we fixed the numerical coefficient 
Ca77 = 1 in the axion-photon coupling for simplicity. Since the axion mass is not given by 
m, we must treat m* and g a77 as different parameters. Cosmological bounds on such models 
are investigated in [61]. The heavy axions which decay at early times lead to various effects 
on the cosmological observations, such as the distortion in the spectrum of CMB, and the 
disagreement of the theoretical prediction with the observed light element abundance. Each 
of the observations gives a severe constraints on the coupling <7 OTr In order to avoid these 
constraints, the lifetime of axions must be short enough to decay earlier than the onset of 
the BBN, i 7 < 10 -2 sec. This requirement leads a bound 

/ F \ ~ 2/3 

S > 4.6 x 1CT 17 x [j^^) ■ (4.38) 

In summary, if the value of S is sufficiently large, the mass of axions is determined by 
the term proportional to S, and this parameter is again constrained by various observational 
considerations. In order to avoid the observational bounds, axions should be light enough 
to avoid the overclosure of the universe, whose condition is given by Eq. (4.34) or (4.36), or 
heavy enough to decay before the BBN epoch, whose condition is given by Eq. (4.38). 



5 Summary and discussions 

In this paper, we have investigated the axion cosmology with iVow > 1) where domain 
walls survive for a long time in the early universe. We performed three-dimensional lattice 
simulations to follow the evolution of string- wall networks in the scaling regime, and computed 
the spectra of axions and gravitational waves produced from them. The radiated axion is 
mildly relativistic with the momentum comparable to the axions mass, which is similar to 
the result of the model with iVow = 1 [11]. The ratio of the mean momentum of axions to 
their mass, given by Eq. (3.8), has a value e a ~ 1. The spectrum of gravitational waves has a 
peak at the momentum comparable to the Hubble scale k/R ~ H, with a plateau extending 
up to the momentum comparable to the inverse of the width of domain walls k/R ~ m. 
The magnitude of gravitational waves is given by Eq. (3.13) with the numerical coefficient 
^gw — 5. 

Using the numerical results, we calculated relic densities of axions and gravitational 
waves produced by long-lived string- wall networks. We found that the significant amounts 
of axions are produced from topological defects, and the requirement that they should not 
overdose the universe gives a stringent constraint on the model parameters. Combined with 
bounds given by astrophysical observations and the experimental result on the NEDM, we 
showed that the most part of parameter regions are excluded unless the value of the phase 
of the bias term S is extremely fine-tuned. In such parameter regions, the amplitude of 
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gravitational waves is too weak to observe even in future gravitational wave detectors with 
high sensitivities. 

Cosmological domain wall problem leads to serious constraints on the axion models. 
The options that we can choose to avoid this problem are listed as follows: 

1. The axion is realized such that ./Vdw = 1- This can be obtained by introducing a 
heavy quark 4 [5]. Domain walls disappear shortly after their formation, and there 
is no cosmological problem. In this case, axion mass is constrained into m ~ 10 -3 - 
l(T 2 eV [11]. 

2. The PQ symmetry is broken during inflation. In this case, the population of topological 
defects is wiped away beyond the horizon scale, and domain walls are not problematic. 
However, the observation of the isocurvature fluctuations in anisotropics of CMB gives 
constraints on the model parameters [15]. 

3. The existence of 3 term annihilates domain walls. This possibility is allowed only if 
the phase 5 of the bias term is extremely suppressed. However, the value of 3 is still 
constrained by cosmology, such as Eqs. (4.34), (4.36), and (4.38). 

If the PQ symmetry is a global symmetry arising accidentally in the low energy effective 
theory, it is easily violated by higher dimensional operators which is suppressed by powers 
of the Planck mass Mp [19]. These operators might be identified as the 3 term, and their 
magnitude must be suppressed in order not to shift the value of 9. The existence of long- 
lived domain walls makes a problem much worse, since the value of 3 cannot be arbitrary 
small. Hence we conclude that the axion models with iVrjw > 1 seems to be excluded from 
cosmological considerations, except for the cases 2 and 3 enumerated above. However, in the 
cases 2 and 3, some tunings on other theoretical parameters such as the inflationary scale 
(for case 2) or the value of 5 (for case 3) are required. Furthermore, even in the case 1, 
there still remains a question of how 3 is suppressed, which should be addressed from some 
fundamental theories. 
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A Numerical schemes 

In this appendix, we describe the analysis method which we used in the numerical studies. 
Our numerical code is the combined version of that we used in past numerical studies [11, 
20, 29, 62]. 

4 This can also be realized by embedding the discrete subgroup Zjv dw into a continuous group [14]. However, 
this requires another continuous (gauge or global) group in addition to the standard model symmetry groups. 



-27- 



A.l Field equations 

The equation of motion for $ is obtained by varying the Lagrangian density given by Eq. (2.1) 
with the potential 

Vm = j(m 2 - V 2 ? + ^ ( 1 - ^ cos N DW e) . (A.l) 
4 jv dw V V J 

The second term of Eq. (A.l) is different from that of Eq. (2.2). The original potential (2.2) 
is not well defined at $ = 0, which causes instabilities in the numerical studies. We avoid this 
singularity by multiplying a factor |$| in front of the cosine term. The difference between 
Eq. (2.2) and Eq. (A.l) is not important in the bulk region on which |$| = rj, and we observe 
that the quantitative behavior such as time evolution of the scaling parameter of topological 
defects is not so much affected by this modification. 

In the numerical studies, we follow the evolution of two real scalar fields (f)\ and fa 
which are defined by <]? = (f>\ + ifo. We normalize all dimensionful quantities in the unit of 
rj, i.e. rj = 1. With this normalization, the equations of motion for two scalar fields in the 
FRW background become 

V 2 

01 + 3F01 



9 

m 



-A0i(0i + 02 - 1) + ^72— (cos 6> cos AW? + A DW 

sin#sinA 7 DW0)> (A-2) 



V 2 

6 2 + 3H<f> 2 



R 2 (ty A 

2 



A<M<Ai + 02 - !) + -r^— (sin9 cos N DW 9 - N DW cos6 smN DW 6). (A.3) 



rn 

If we use the conformal time r (dr = dt/R(t)), the above equations become 

/l'-V 2 /! 

-A/i(|/i| 2 + |/ 2 | 2 - -R 2 ) + ^-^(costfcos N DW 9 + A DW sin^sinA D w^), (A.4) 



fl -v 2 / 2 



,2 

-A/ 2 (|/i| 2 + I/2I 2 - # 2 ) + — (sin cos Adw$ - Afow cos sin , (A.5) 

iV DW 

where / = ii^ = f 1+1/2, and a prime denotes a derivative with respect to the conformal time 
r. Here, we assume the radiation dominated background. This model has three parameters: 
A, Adw and m. In the numerical studies, we fix the values of A and m, and vary the value 
of AW- 

A. 2 Discretization of differential equations 

To execute the time evolution which follows from Eqs. (A.4) and (A.5), we used the symplectic 
integration scheme developed in [63]. In this numerical scheme the (implicit) Hamiltonian of 
the system is exactly conserved and there is no secular growth of the conserved quantity. In 
practice, we compute the following sequence of mappings, 

/(r, x) -> f(r + c r At, x) = /(r, x) + c r Arf'(T, x), (A.6) 
f'(r, x) -> f'(r + d r Ar, x) = f'(r, x) + d r Ar/"(r, x), (A.7) 
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from r = 1 to r = 4, in each time step. The values of coefficients c r and d r are given by [63] 



1 _ 1 - 2 1 / 3 

Cl ~ Q ~ 2(2 - 2V3) ' C2 - c 3 - 2 ( 2 _ 2 i/3) ' 

1 2 x /3 

It turns out that this integration scheme has fourth-order accuracy in time. 

The spatial derivative and the Laplacian of the scalar fields are computed by using the 
finite difference method. For example, the Laplacian is given by 



(V 2 'f)i,j,k — 10 ^ A x 2 [16(/i+ij,fe + fi-l,j,k + fi,j+l,k + fi,j-l,k + fi,j,k+l + fi,j,k-l 



1 

12(Ax) 

(fi+2,j,k + fi-2,j,k + fi,j+2,k + fi,j-2,k + fi,j,k+2 + fi,j,k-2) 

90/ij,fc] , (A.9) 



where /i j jfc = /(r, Xj, y.,-, z^), and Ax = L/iV is the grid spacing. This formula has fourth- 
order accuracy in Ax. We imposed the periodic boundary condition on the boundaries of 
the simulation box. 



A. 3 Initial conditions 

We generate initial conditions in the similar way with Refs. [20, 35]. We treat <j)\ and 4>2 
as two independent real scalar fields, and assume that each of them has massless quantum 
fluctuations at the initial time. With this assumption, we generate initial conditions as 
Gaussian random amplitudes in momentum space. Then we Fourier transform them into 
the configuration space to give the spatial distribution of the fields. We emphasize that this 
choice of initial conditions does not take account of the correct circumstances in the QCD 
epoch. However, we expect that the result in the late time of the simulations are qualitatively 
unchanged if we used the different initial conditions, since the evolution of defect networks 
in the scaling regime is not so much affected by the initial field configurations. See appendix 
of Ref. [20] for more discussions. 



B Identification of topological defects 

Using the lattice algorithm described in Appendix A, we can solve the time evolution of two 
real scalar fields </>i(i, x) and (p2(t,x) in the comoving simulation box. From these data of 
scalar fields, we estimate various physical quantities, such as the length of strings, the area of 
domain walls, and the spectra of axions and gravitational waves radiated by defect networks. 
In order to calculate such physical quantities, we have to identify the position of topological 
defects in the simulation box. In this appendix, we summarize our identification method of 
topological defects, and describe how to calculate scaling parameters of strings £ and domain 
walls A defined in Sec. 3.1. 

B.l Identification of strings 

We used the method developed in [29] to identify the position of strings. Consider a quadrate 
which consists of four neighboring grids in the simulation box (see points A, B, C, and D in 
Fig. 12). At each vertex point, we estimate the phase of the complex scalar field 6 from the 
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data of two real scalar fields 4>\ and 02- Then, we obtain a mapping from the quadrate in the 
real space to the phase distribution in the field space, as shown in Fig. 12. Let us denote the 
minimum range of the phase which contains all the images of four vertices as A9. The string 
exists within the quadrate if A9 > ir is satisfied, as long as the value of the phase changes 
continuously around the core of strings. 

We can also determine the position at which a sting penetrates the quadrate. At the 
boundary of the quadrate, we determine the positions of points on which <pi = or = 
is satisfied, by using linear interpolation of the values of (f>i and fo- If the string correctly 
penetrate the quadrate, there are two points with (j>\ = and two points with = on 
the boundary of the quadrate. Then we determine the position of the string as the point 
at which the line connecting two points with <f>\ = intersects with the line connecting two 
points with 02 = (see left-hand panel of Fig. 12). 

This simple criterion breaks down when the quadrate is penetrated by more than two 
strings. However, we observed that such region is at most 1% of the whole simulation 
box when the system relaxes into the scaling regime. Therefore, this scheme enables us to 
determine the position of strings with at least 99% accuracy. 

B.2 Identification of domain walls 

The identification of domain walls is simpler than that of strings. We separate the region 
of the phase of < 9 < 2ir, into iV^w domains. For example, if JVdw = 3, we obtain 
three domains with < 9 < tt/3 and 5ir/3 < 9 < 2ir (vac. 1), tt/3 < 9 < ir (vac. 2), and 
7r < 9 < 5ir/3 (vac. 3), as shown in Fig. 12. At each grid point, we compute the phase of 
the scalar field and assign the number of the vacuum domain (i.e. vac. 1, vac. 2, or vac. 3) 
which contains that point. Let us call the neighboring grid points that differ by a unit lattice 
spacing Ax = L /N as the "link" . We identify that the domain wall intersects the link if the 
number of vacuum is different between the two ends of the link. 

B.3 Calculation of scaling parameters 

The scaling parameters of strings and domain walls are given by 

e= A*ring t 2 and A= P™U tf (B.l) 
Mstring ^wall 

where /String is the energy density of strings, p wa u is the energy density of domain walls, 
Mstring is the mass energy of strings per unit length, and o" wa n is the surface mass density 
of domain walls. It might be difficult to determine the values of String and <r wa n directly 
from the data of scalar fields since they are quantities obtained, by integrating /String 

and 

Pwaii ° ver the coordinates which are perpendicular to the direction along which the defects 
extend. Instead of estimating /i s t r ing and <r wa u, we calculate the length / of strings and area 
A of domain walls in the comoving coordinates, such that 

/^string 

Pstnng-^7 and PwaU = ^F, (B.2) 

where V = L 3 is the comoving volume of the simulation box. Then, the scaling parameters 
are given by 

It 2 At 
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Figure 12. Schematics of the identification method of topological defects. Left panel shows a 
quadrate with four vertex points (A, B, C, and D) penetrated by a string and domain walls in the 
real space. Right panel shows a mapping of them in the field space. In this figure, we assume the 
model with TVdw = 3. Hence, there are three domain walls attached to the string. The blue (wall 1), 
pink (wall 2), and yellow (wall 3) lines represent the locations of the center of domain walls, which 
correspond to 9 — 7r/3, 9 = n, and 9 = 57r/3, respectively, in the field space. The region is separated 
into three domains, vac. 1, vac. 2, and vac. 3, which are surrounded by domain walls on the boundaries. 
The dashed lines correspond to the loci of 4>i = and (j>2 = 0, which intersect on the core of the 
string. In the field space, we define A9 as the minimum range of the phase which contains all the 
images of four vertices. A string penetrates the quadrate if the condition A9 > ir is satisfied. 



The length of strings I can be computed by connecting the loci of strings identified by 
the method described in the previous subsection. The computation of the area of domain 
walls A is not so straightforward, since it depends on the way how we define the segment 
of A in each lattice. Here, we use the simple algorithm introduced by Press, Ryden, and 
Spergel [64]. Define the quantity 6± which takes the value 1 if the number of vacuum is 
different between the two ends of the link, and the value otherwise. We sum up 5± for each 
of the links with the weighting factor 

where AA = (Ax) 2 is the area of one grid surface, and 6 x , etc. is a derivative of #(x) with 
respect to x. The weighting factor given by the function of V# gives the average number of 
links per area segment. This method takes account of the orientation of the surface of walls 
in the lattice. See [64] for details. 

C Calculation of power spectrum of radiated axions 

The calculation of power spectrum of axions radiated from topological defects is executed by 
using the masking analysis method developed in Refs. [11, 29]. In this appendix, we describe 
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the main line of the calculation scheme. Readers may also refer to [11, 29] for complete 
description on the analysis method. 

The power spectrum P(k, t) of axions is defined by 

l(d(i,k)*d(t,k')) = ^«5(3)( k -k')P(M), (CI) 

where (...) represents an ensemble average. a(t, k) is the Fourier component of the time 
derivative of the axion field 

d(t,k) = J d 3 xe ik x d(t,x), (C.2) 
where the value of d(t, x) is obtained from the simulated data of $ and $ 



a(t, x) = Im 



(C.3) 



If a(t,x) is a free field, it can be shown that P(k,t) is proportional to the energy spectrum 
of axions [11]. 

The data of a obtained by numerical simulations contain the contaminations from the 
core of strings or domain walls. We mask the contribution of axion field near strings or 
domain walls by introducing a window function 

TI ,/ s f (near defects) ,< 
W & = { 1 (elsewhere) " (C ' 4) 

Then, we compute the power spectrum P(k) of the masked axion field 

d(x) = W(x)d(x). (C.5) 

Unfortunately, the masked spectrum P(k) is not equivalent to the true spectrum of free 
axions. This discrepancy can be resolved by using the pseudo-power spectrum estimator 
(PPSE) 

k 2 r rlk' 

PppsE(k) = - J ^M~\k, k')P(k% (C6) 
where M~ 1 (k, k') is a matrix which satisfies 

r k' 2 rlk' 9tt 2 

J ^T M " 1(Jfe ' k')M(k', k") = ^b(k - k"). (C.7) 
Here, M(k, k') is the window weight matrix defined by 

where W(k) is the Fourier transform of the window function W(x), Vtj, is a unit vector 
representing the direction of k, and d£lk = d cos 9d(j). It turns out that the ensemble average 
of PppsE(k) is equivalent to the spectrum of radiated axions [29]. 

Naively, M{k,k') in Eq. (C.8) can be evaluated by averaging |W(k — k')| 2 /l/ 2 over 
discrete k and k' in the lattice. This however requires 0(N 2 ) arithmetics, which is practically 
prohibited when the grid number N is large as is in our analysis. Instead, we approximate 
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M(k, k') by partially taking the continuous limit, which we here describe briefly. With a 
transformation P = k — k' and R = k t> k , Eq. (C.8) can be rewritten as 



{^) 2 V 2 k 2 k U 



d 3 P\W{P)\' 



M(k,k') 



In the continuous limit, the integration over R can be performed analytically. 
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where P = |P|, R= |R| and \i = ^p-- We thus obtain 
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Evaluation of Eq. (C.ll) in a lattice has a good scaling with only O(N) arithmetics. On the 
other hand, the approximation can cause a systematic error in estimation of -Pppse(^) 5 which 
shall be discussed shortly later. 

The power spectrum PppsE(k,t) contains the contributions of axions produced before 
the time t. This means that it contains the contributions of axions radiated at initial epoch 
where the defects do not relax into the scaling regime. Since we are interested in the spec- 
trum of axions radiated from topological defects which evolve along the scaling solution, we 
must subtract the contribution of axions produced at initial relaxation regime. This can 
be achieved by evaluating the power spectra for selected time steps t\ and t 2 > t±, and 
subtracting the spectrum obtained at £2 from that obtained at ti [11] 

3 



AP(k, t 2 ) = Ppp S E(k, t 2 ) - PppsE{k, ti) 



Lo{k,t 2 ) fR(h 



(C.12) 



u(k,h) \R(t 2 ) 

where (io(k,t 2 )/uj(k,ti))(R(ti)/R(t 2 ))' i is a dilution factor which comes from the effect of the 
cosmic expansion, and 

u(k,t) = ^m 2 + k 2 /R{t) 2 (C.13) 

is the energy of axions with momentum k/R(t). 

Finally, we comment on the estimation of errors of the power spectrum. Error is esti- 
mated by the following procedure based on the Monte Carlo sampling. After we compute 
PppsE^k), we generate d(k) as Gaussian random field with a variance PppsE(k)/k 2 in the 
Fourier space. Then we execute the inverse Fourier transformation of d(k) to obtain the data 
of d(x) in the real space. Let us denote this spurious field data as dMc( x )- We apply the 
same masking procedure to the field dMc( x ) and obtain the power spectrum -Pppse,mc(&)- 
We iterate this procedure by hmc times, and compute the mean power spectrum 

^ "MC 

PPPSE,MC(A:) = P PPSE,MC(*0 (C14) 
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and the covariance matrix of -Pppse,mc(^) 



-, HMO 

C(k, k') = ^(^P S E,MC(^) " ^PPSE,Mc(A ; ))(P^psE,Mc(fc / ) " JW.MC (&')), (C15) 

where -fppsEMc(^) ^ s ^ ne P ower spectrum obtained at i-th sampling step. Although PPSE 
should be unbiased as formally shown in Ref. [29], we observed that -Pppse,mc(&) does not 
converge into -Pppse(^) obtained by the field-theoretic simulation. This bias (-Pppse,mc(&) — 
-Pppse(^)) arises mainly from use of the approximation M(k, k') in Eq. (C.ll). In our analysis, 
we regard the bias as the systematic error. Then the covariance matrix C(k,k'), which 
accounts for both the statistical and systematic errors, can be given by 

C(k, k') = C(k, k') + (Pppse.mcW - PppseWX^ppse.mc^O - PppsE(k')) 

= £(^PSE,Mc(*0 " J PpPSE(fc))(P^p S E,Mc(^) " ^PPSE^'))- (C16) 

«MC 

4 = 1 

In addition to these errors in estimation of -Pppse(^) in each realization, -Pppse(^) also varies 
from one realization to another, variance from which should be also included as a statistical 
error. Thus energy spectrum averaged over realizations and its covariance matrix can be 
finally given by 



fppsE(fc) = ^E P ppseW> (C.17) 

T ' 3=1 

C(k,k') = -Lj2 GU) (k,k') 

1 

+^3T D P ppse(*0 - PppsE(k)][P^ SE (k') - P PPSE (k% (C.18) 



3=1 



N r 

3=1 

where N r is the number of realizations, and the subscript (J) indicates a realization. We 
execute the numerical simulations with N r = 20 different realizations of initial field configu- 
rations. The results we present in the text is obtained from -Pppse(^) an d C(k,k'). 

D Calculation of gravitational waves 

We use the method introduced in [65] to calculate the spectrum of gravitational waves. The 
evolution of gravitational waves is described by spatial metric perturbations hij around the 
FRW background 

ds 2 = -dt 2 + R 2 {t)(5ij + h i j)da*da i . (D.l) 
Let us define the proper fluctuations as 

Xi i (x,r) = i?(r)/ iii (x,r), (D.2) 

and denote their Fourier components as 

Xy(k,r) = J d^^Xiji^r). (D.3) 
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In the Fourier space, the equations of motion for metric perturbations reduce to the form [35] 



dy 2 



+ 1 



Xij 



S„ 



k 2 ' 



where we defined the variable y = kr and the source 



Sij(k,T) = 167rG J R(r)^ T (k,r). 



(D.4) 



(D.5) 



The transverse traceless (TT) part of the stress-energy tensor is computed by applying the 
projection operator 



rg T (k,r)=A ii)fc ,(fc)r -(k,r) 

= A ijM (k){d k ^*d^}(k,r), 

A ij>k i(k) = P ik (k)Pji(k) - ^Pij(k)P kl (k), 



Pijik) — Sij k{kj , 



(D.6) 
(D.7) 
(D.8) 



where k = k/|k|, and {d k <&*<9/<I>} (k, r) is the Fourier transform of <9fc<£*(x, r)c^<i>(x, r). 

There are two homogeneous solutions of Eq. (D.4), cosy and siny. The special solution 
is given by the time integral of the source term convoluted with two homogeneous solutions 



Xij (k, r) = C$ (k, r) cos kr + c\f (k, r) sin kr, 



where 



and 



v , . , S i3 -(k,r') _ 16ttG- (1) 
d y smfer = - i[ 5-C i V(k,r) I 



A- 2 



A- 2 



dy' cos kr , 9 = , 9 C^ ; (k,r). 



A- 2 



A; 2 



The energy density of the gravitational waves is given by 

1 



Pgw(r) 



32itGR 4 ^ 



(Xiji*, 7-)x'y(x, r)) 



(D.9) 
(D.10) 
(D.ll) 

(D.12) 



where a prime denotes a derivative with respect to conformal time r. We replace the ensemble 
average in Eq. (D.12) by an average over a volume V of the comoving simulation box, 



Pgw(r) 



1 



1 



d 3 k 



32ttGR 4 V J (2vr) 



X / «(k,r)xS(k,r). 



\3 A y 



(D.13) 



Substituting the solution given by Eq. (D.9), and ignoring the terms with higher order in 
RH, we obtain 



4ttG f d 3 k 1 



R*V J (2tt) 3 A: 2 



£ I*. 



+ 



C 



(2) 



(D.14) 



where we averaged over a period of the oscillation of x?j(k, r) with time. The spectrum of 
gravitational waves is given by 



dp, 



gw 



Gk 



dink 2tt 2 VR 4 



/^e(i^ 3 



^(2) 



(D.15) 
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We define the dimensionless energy spectrum of gravitational waves 



dp gw /d In k 
Pc(r) 



G 2 



3irV R A H 2 



S k (r) 



(D.16) 



where 



s k (r) = k J dntY^ncW 



+ 



c 



(2) 



(D.17) 



and Pc(t) is the critical density of the universe at conformal time r 



Pc(r) 



3H 2 



(D.18) 



In radiation dominated background, we expect that f2 gw (/c,T) oc S k {r), since R A H 2 remains 
constant. Therefore, we calculate Sk{r) in the numerical studies. Note that f2 gw (/c,r) given 
by Eq. (D.16) does not correspond to the spectrum observed today. The spectrum of grav- 
itational waves observed today is obtained by multiplying Q gw (fc,r) with a damping factor 
which arises from the following matter-dominated epoch. 

As we discussed in Sec. 3.3, we subtract the contribution of gravitational waves produced 
at the initial epoch where the defects do not relax into the scaling regime. This subtraction 
can be executed as follows. Since the energy density of gravitational waves shifts as oc i?~ 4 , 
the quantity ^(r) oc R(t) dp gw (T) / din k is proportional to the number of gravitons with 
the momentum k per comoving volume. Therefore, the difference of the spectral function 



AS k (r) = S fc (T) - S k (r g ) 



(D.19) 



gives the spectrum of gravitational waves produced during the time interval [t s ,t]. Here T g 
is a reference time which we treat as a free parameter in the numerical simulations. 

E Surface mass density of domain walls 

In this appendix, following Ref. [8], we shortly comment on the surface mass density of domain 
walls in the models used in numerical simulations. First, consider the model of scalar field (j> 
with the potential given by 



v{4>) 



T 1 



-2\2 



Vr 



(E.l) 



where A r and rj r are constants. This model has a discrete Zi symmetry in which the scalar 
field transforms as <j> — > —(f). In the flat Minkowski background, the classical field equation 
of this model has a static solution 



rj r tanh 



rj r z 



(E.2) 



which describes a domain wall lying perpendicular to the z-axis. By integrating the energy 
density of the field over the direction perpendicular to the surface of the wall, we obtain the 
surface mass density 
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Next, let us consider the axionic model with the Lagrangian given by Eq. (2.1). If we 
restrict ourselves to low energy configurations at the QCD scale, we can put = rje 18 , which 
gives the effective potential for 6 



2 2 2 

rf ,„ > 2 m V 



A> = + -2^(1- cos iV DW 0). (E.4) 

1 JV DW 

The classical field equation in the Minkowski background again yields a planner solution 
perpendicular to the z-axis 

°( z ) = ^J— + tan _1 exp(mz), k = 0, 1, . . . , iV DW - 1. (E.5) 

Integrating out the energy density, we obtain the surface mass density of the domain wall 

Comparing Eqs. (E.3) and (E.6), the values of cr wa ii coincide between these two models 
if we set A r = 4m 2 /3rj 2 , r] r = (3/2) 1//2 ry, and Afow = 2. Also, the width of domain walls 
becomes the same order 6 W ~ {\ T /2)~ l / 2 ri~ 1 ~ mT 1 . Substituting these values of A r and ij r 
into Eq. (E.l), we obtain the potential given by (3.3). 

We note that the surface mass density of axionic domain wall, given by Eq. (E.6) should 
be modified if we include the structure of the neutral pion field which varies inside the wall. 
In this case, the surface mass density is given by [66] 

cj wall = 4.32/^771^ /N DW ~ 9.23m V 2 /N 2 w , (E.7) 

where f n is the pion decay constant, and m n is the mass of the pion. In the second equality, 
we used rriu/md ~ 0.48 [3], where m u and are the mass of up quark and down quark, 
respectively. 
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